Antje Hofmann: Learning Syntactic Relations with the Mole Task

RePsychLing in SMLP2022

Author

Reinhold Kliegl

Published

August 24, 2022

1 Background

Children’s (age: 4-8 years)reaction times in a task teaching them syntactic relations.

1.1 Overview

  • Original analysis is by Antje Hofmann.
  • MixedModels.jl version
  • Addition of new chunks illustrate
    • selection of parsimonious LMM using random-effects PCA
    • plotting conditional means
    • illustration of borrowing strength

2 Readme

2.1 Design

2.2 Variables

  • Subj: Participant ID (renamed from ID; random factor)
  • Item: Word ID (random factor)
  • age: 4 - 8 years
  • Block (within-Subj/within-Item):
    • 1st Learning
    • 2nd Learning
    • Disruption
    • Recovery
  • Target(renamend fom targetness)
    • non-syllable target
    • syllable target
  • rt: response time

3 Setup

3.1 Packages

First attach the MixedModels.jl package and other packages for plotting. The CairoMakie.jl package allows the Makie graphics system [@Danisch2021] to generate high quality static images. Activate that package with the SVG (Scalable Vector Graphics) backend.

Code
using AlgebraOfGraphics
using Arrow
using CairoMakie       # graphics back-end
using CategoricalArrays
using Chain
using DataFrames
using DataFrameMacros  # simplified dplyr-like data wrangling
using MixedModels
using MixedModelsMakie # diagnostic plots
using ProgressMeter
using Random           # random number generators
using RCall            # call R from Julia
using StatsModels

using AlgebraOfGraphics: boxplot
using AlgebraOfGraphics: density

using MixedModelsMakie: qqnorm
using MixedModelsMakie: ridgeplot
using MixedModelsMakie: scatter
using MixedModelsMakie: caterpillar

ProgressMeter.ijulia_behavior(:clear);
CairoMakie.activate!(; type="svg");
  • The data are available as an arrow file.
  • Most of preprocessing was done with R in RStudio (see Hofmann_Maulwurf.Rmd).
  • Order of factor levels should be checked.
Code
dat = DataFrame(Arrow.Table("./data/Hofmann_Maulwurf_rt.arrow"))
transform!(dat,
    :Target => categorical => :Target,
    :Block  => categorical => :Block,
    :age => (x -> x .- 6) => :a1, # center age (linear) at 6 years
    :rt => (x -> log.(x)) => :lrt)
describe(dat)

8 rows × 7 columns

variablemeanminmedianmaxnmissingeltype
SymbolUnion…AnyUnion…AnyInt64Type
1Subjp01p850Union{Missing, String}
2Item01PA.ogg32SG.ogg0Union{Missing, String}
3age6.210813.76.318.040Union{Missing, Float64}
4Block1st LearningRecovery0Union{Missing, CategoricalValue{String, UInt32}}
5TargetNon-target syllableTarget syllable0Union{Missing, CategoricalValue{String, UInt32}}
6rt2335.39228.02100.08389.00Union{Missing, Float64}
7a10.210813-2.30.312.040Float64
8lrt7.673675.429357.649699.034680Float64
  • Centering age at six years yields an interpretable GM
  • Factor levels can also be set when contrasts are defined (see below).
  • BoxCox check showed that reaction time rt [ms] should be transformed to speed [1/s] = [Hz]
  • Indicator variables for Target and Block generated in R.

4 LMM analysis

4.1 Contrasts

Code
contrasts = merge(
      Dict(:Target => EffectsCoding()),
      Dict(:Block => SeqDiffCoding()),
      Dict(nm => Grouping() for nm in (:Subj, :Item))
   );

4.2 Varying only intercepts LMM m_voi

Code
f_voi1 = @formula(lrt ~  1 + Block * Target * a1 + (1 | Subj) + (1 | Item));
m_voi1 = fit(MixedModel, f_voi1, dat; contrasts)
Minimizing 65    Time: 0:00:01 (17.57 ms/it)
Est. SE z p σ_Item σ_Subj
(Intercept) 7.6895 0.0371 207.05 <1e-99 0.0239 0.2637
Block: 2nd Learning -0.0920 0.0104 -8.87 <1e-18
Block: Disruption 0.0516 0.0126 4.09 <1e-04
Block: Recovery -0.0477 0.0125 -3.81 0.0001
Target: Target syllable -0.0108 0.0040 -2.69 0.0072
a1 -0.0896 0.0324 -2.77 0.0056
Block: 2nd Learning & Target: Target syllable -0.0057 0.0104 -0.55 0.5854
Block: Disruption & Target: Target syllable -0.0262 0.0124 -2.11 0.0347
Block: Recovery & Target: Target syllable 0.0020 0.0123 0.16 0.8741
Block: 2nd Learning & a1 0.0347 0.0091 3.83 0.0001
Block: Disruption & a1 0.0214 0.0108 1.97 0.0486
Block: Recovery & a1 -0.0159 0.0108 -1.47 0.1409
Target: Target syllable & a1 -0.0052 0.0035 -1.46 0.1430
Block: 2nd Learning & Target: Target syllable & a1 -0.0113 0.0090 -1.25 0.2111
Block: Disruption & Target: Target syllable & a1 0.0113 0.0108 1.05 0.2933
Block: Recovery & Target: Target syllable & a1 0.0086 0.0107 0.80 0.4248
Residual 0.3004

4.3 Extract indicator variables

Code
X = modelmatrix(m_voi1)
dat.trng = X[:,2];
dat.drpt = X[:,3];
dat.rcvr = X[:,4];
dat.trgt = X[:,5];
describe(dat)

12 rows × 7 columns

variablemeanminmedianmaxnmissingeltype
SymbolUnion…AnyUnion…AnyInt64Type
1Subjp01p850Union{Missing, String}
2Item01PA.ogg32SG.ogg0Union{Missing, String}
3age6.210813.76.318.040Union{Missing, Float64}
4Block1st LearningRecovery0Union{Missing, CategoricalValue{String, UInt32}}
5TargetNon-target syllableTarget syllable0Union{Missing, CategoricalValue{String, UInt32}}
6rt2335.39228.02100.08389.00Union{Missing, Float64}
7a10.210813-2.30.312.040Float64
8lrt7.673675.429357.649699.034680Float64
9trng-0.014562-0.750.250.250Float64
10drpt-0.0528019-0.5-0.50.50Float64
11rcvr0.0524467-0.25-0.250.750Float64
12trgt0.00331492-1.01.01.00Float64

Switch to indicator variables and refit the model.

Code
f_voi2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1 + (1 | Subj) + (1 | Item));
m_voi2 = fit(MixedModel, f_voi2, dat; contrasts)
Est. SE z p σ_Item σ_Subj
(Intercept) 7.6895 0.0371 207.05 <1e-99 0.0239 0.2637
trng -0.0920 0.0104 -8.87 <1e-18
drpt 0.0516 0.0126 4.09 <1e-04
rcvr -0.0477 0.0125 -3.81 0.0001
trgt -0.0108 0.0040 -2.69 0.0072
a1 -0.0896 0.0324 -2.77 0.0056
trng & trgt -0.0057 0.0104 -0.55 0.5854
drpt & trgt -0.0262 0.0124 -2.11 0.0347
rcvr & trgt 0.0020 0.0123 0.16 0.8741
trng & a1 0.0347 0.0091 3.83 0.0001
drpt & a1 0.0214 0.0108 1.97 0.0486
rcvr & a1 -0.0159 0.0108 -1.47 0.1409
trgt & a1 -0.0052 0.0035 -1.46 0.1430
trng & trgt & a1 -0.0113 0.0090 -1.25 0.2111
drpt & trgt & a1 0.0113 0.0108 1.05 0.2933
rcvr & trgt & a1 0.0086 0.0107 0.80 0.4248
Residual 0.3004

They are equivalent.

4.4 A zero-correlation parameter LMM m_zcp

Code
f_zcp1 = @formula(lrt ~  1 +  Block * Target * a1 + zerocorr(1 + Block * Target | Subj) + (1 + a1 | Item));
m_zcp1 = fit(MixedModel, f_zcp1, dat; contrasts);

show(issingular(m_zcp1))
VarCorr(m_zcp1)
Minimizing 463   Time: 0:00:00 ( 1.23 ms/it)
false
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0742252 0.2724431
Block: 2nd Learning 0.0362305 0.1903431 .
Block: Disruption 0.0156482 0.1250926 . .
Block: Recovery 0.0133091 0.1153650 . . .
Target: Target syllable 0.0002450 0.0156523 . . . .
Block: 2nd Learning & Target: Target syllable 0.0001921 0.0138594 . . . . .
Block: Disruption & Target: Target syllable 0.0032728 0.0572082 . . . . . .
Block: Recovery & Target: Target syllable 0.0006683 0.0258522 . . . . . . .
Item (Intercept) 0.0009504 0.0308293
a1 0.0006760 0.0259993 -0.83
Residual 0.0767202 0.2769841

Again, check the equivalence.

Code
f_zcp2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1  +
                zerocorr(1 + (trng+drpt+rcvr) * trgt | Subj) + (1 + a1 | Item));
m_zcp2 = fit(MixedModel, f_zcp2, dat; contrasts);

show(issingular(m_zcp2))
VarCorr(m_zcp2)
Minimizing 463   Time: 0:00:00 ( 0.53 ms/it)
false
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0742252 0.2724431
trng 0.0362305 0.1903431 .
drpt 0.0156482 0.1250926 . .
rcvr 0.0133091 0.1153650 . . .
trgt 0.0002450 0.0156523 . . . .
trng & trgt 0.0001921 0.0138594 . . . . .
drpt & trgt 0.0032728 0.0572082 . . . . . .
rcvr & trgt 0.0006683 0.0258522 . . . . . . .
Item (Intercept) 0.0009504 0.0308293
a1 0.0006760 0.0259993 -0.83
Residual 0.0767202 0.2769841

4.5 A complex parameter LMM m_cpx

Code
m_cpx = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + trgt * (trng+drpt+rcvr) | Subj) + (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_cpx))    # not ok
show(m_cpx.PCA[:Subj])     # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_cpx))
VarCorr(m_cpx)
Minimizing 1355      Time: 0:00:00 ( 0.64 ms/it)
true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .      .
 trgt         -1.0    1.0     .      .      .      .      .      .
 trng          0.35  -0.35   1.0     .      .      .      .      .
 drpt          0.06  -0.06   0.35   1.0     .      .      .      .
 rcvr         -0.42   0.42  -0.23  -0.48   1.0     .      .      .
 trgt & trng   0.26  -0.26  -0.12   0.09  -0.16   1.0     .      .
 trgt & drpt  -0.4    0.4   -0.07   0.14  -0.01  -0.76   1.0     .
 trgt & rcvr   0.61  -0.61   0.16   0.04  -0.3    0.48  -0.85   1.0

Normalized cumulative variances:
[0.441, 0.667, 0.8072, 0.906, 0.9581, 1.0, 1.0, 1.0]

Component loadings
 
               PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8
 (Intercept)  -0.46   0.14   0.36  -0.2   -0.28   0.12   0.05  -0.71
 trgt          0.46  -0.14  -0.36   0.2    0.28  -0.12  -0.05  -0.71
 trng         -0.19   0.43   0.17   0.76  -0.0   -0.4   -0.13   0.0
 drpt         -0.09   0.49  -0.59   0.15  -0.19   0.57   0.13  -0.0
 rcvr          0.26  -0.42   0.31   0.47  -0.39   0.51  -0.15   0.0
 trgt & trng  -0.31  -0.36  -0.49  -0.01  -0.52  -0.34  -0.38   0.0
 trgt & drpt   0.39   0.44   0.14  -0.29  -0.16   0.02  -0.72   0.0
 trgt & rcvr  -0.46  -0.18  -0.04   0.08   0.6    0.33  -0.53   0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2393.5710                         
[2]         56  2353.2482  40.3229      28  0.0619
──────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0757920 0.2753034
trgt 0.0001016 0.0100800 -1.00
trng 0.0334990 0.1830274 +0.35 -0.35
drpt 0.0198176 0.1407751 +0.06 -0.06 +0.35
rcvr 0.0191682 0.1384494 -0.42 +0.42 -0.23 -0.48
trgt & trng 0.0020614 0.0454021 +0.26 -0.26 -0.12 +0.09 -0.16
trgt & drpt 0.0145492 0.1206203 -0.40 +0.40 -0.07 +0.14 -0.01 -0.76
trgt & rcvr 0.0085469 0.0924494 +0.61 -0.61 +0.16 +0.04 -0.30 +0.48 -0.85
Item (Intercept) 0.0000000 0.0000000
a1 0.0004752 0.0217991 +NaN
Residual 0.0768486 0.2772159

The deviance improves, but we end up with an overparameterized LMM. Remove the two small VCs of the three Block x Target interactions (i.e., the third contrast of the Block factor).

4.6 A parsimonious parameter LMM m_prm

We remove one of the VC for trgt * trng contrast interaction, that is one of the three interaction terms.

Code
m_prm1 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + trgt *(drpt + rcvr) + trng  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm1))  # not ok
show(m_prm1.PCA[:Subj])   # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm1, m_cpx))
VarCorr(m_prm1)
Minimizing 1465      Time: 0:00:00 ( 0.61 ms/it)
true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .
 trgt         -0.58   1.0     .      .      .      .      .
 drpt          0.04  -0.23   1.0     .      .      .      .
 rcvr         -0.41   0.21  -0.47   1.0     .      .      .
 trgt & drpt  -0.42   0.3    0.21  -0.06   1.0     .      .
 trgt & rcvr   0.65  -0.22   0.03  -0.29  -0.87   1.0     .
 trng          0.33  -0.32   0.32  -0.21  -0.12   0.22   1.0

Normalized cumulative variances:
[0.4076, 0.6504, 0.7761, 0.8828, 0.9634, 1.0, 1.0]

Component loadings
                PC1    PC2    PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.5    0.02  -0.05  -0.36  -0.38  -0.62  -0.3
 trgt          0.38   0.14   0.56   0.49  -0.35  -0.3   -0.26
 drpt         -0.15  -0.61   0.15   0.25   0.61  -0.36  -0.14
 rcvr          0.29   0.44  -0.56   0.23   0.27  -0.52   0.1
 trgt & drpt   0.4   -0.49  -0.05  -0.23  -0.33  -0.3    0.59
 trgt & rcvr  -0.49   0.29   0.35   0.24   0.07  -0.13   0.68
 trng         -0.3   -0.29  -0.47   0.64  -0.42   0.15  -0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + drpt + rcvr + trgt & drpt + trgt & rcvr + trng | Subj) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.5710                          
[2]         48  2334.2053   59.3657      20  <1e-05
[3]         56  2353.2482  -19.0428       8     NaN
───────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0759692 0.2756251
trgt 0.0003130 0.0176929 -0.58
drpt 0.0194377 0.1394192 +0.04 -0.23
rcvr 0.0188254 0.1372057 -0.41 +0.21 -0.47
trgt & drpt 0.0109856 0.1048121 -0.42 +0.30 +0.21 -0.06
trgt & rcvr 0.0075846 0.0870893 +0.65 -0.22 +0.03 -0.29 -0.87
trng 0.0354738 0.1883450 +0.33 -0.32 +0.32 -0.21 -0.12 +0.22
Item (Intercept) 0.0008077 0.0284198
a1 0.0006179 0.0248571 -0.94
Residual 0.0761744 0.2759972

We don’t lose goodness of fit, but are still overparameterized. We remove CPs for Target.

Code
m_prm2 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) +
                  zerocorr(0 + trgt  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm2))  #  ok
show(m_prm2.PCA[:Subj])   #  ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm2,  m_cpx))
VarCorr(m_prm2)
Minimizing 740   Time: 0:00:00 ( 0.55 ms/it)
false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt          0.04   1.0     .      .      .      .     .
 rcvr         -0.41  -0.47   1.0     .      .      .     .
 trng          0.33   0.34  -0.22   1.0     .      .     .
 drpt & trgt  -0.36   0.26  -0.1   -0.07   1.0     .     .
 rcvr & trgt   0.59  -0.0   -0.27   0.12  -0.87   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.3444, 0.5975, 0.7404, 0.8569, 0.9478, 0.9939, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.52   0.07  0.0   -0.04   0.61  -0.58   0.15
 drpt         -0.09   0.61  0.0   -0.05  -0.62  -0.47   0.11
 rcvr          0.3   -0.49  0.0    0.52  -0.16  -0.58  -0.17
 trng         -0.26   0.36  0.0    0.83   0.09   0.3   -0.08
 drpt & trgt   0.47   0.45  0.0   -0.1    0.37  -0.14  -0.65
 rcvr & trgt  -0.59  -0.22  0.0   -0.14  -0.27   0.03  -0.72
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.5710                          
[2]         42  2341.3452   52.2258      14  <1e-05
[3]         56  2353.2482  -11.9030      14     NaN
───────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0760175 0.2757128
drpt 0.0193591 0.1391371 +0.04
rcvr 0.0188566 0.1373195 -0.41 -0.47
trng 0.0351279 0.1874243 +0.33 +0.34 -0.22
drpt & trgt 0.0099356 0.0996773 -0.36 +0.26 -0.10 -0.07
rcvr & trgt 0.0066342 0.0814504 +0.59 -0.00 -0.27 +0.12 -0.87
trgt 0.0002660 0.0163107 . . . . . .
Item (Intercept) 0.0008269 0.0287560
a1 0.0006738 0.0259576 -0.93
Residual 0.0762049 0.2760523

Check the Item-related CP. It is very large, might be spurious.

Code
m_prm3 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  +
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) +
                  zerocorr(0 + trgt  | Subj) +
                  zerocorr(1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm3))  # ok
show(m_prm3.PCA[:Subj])   # ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm3, m_prm2, m_cpx))
VarCorr(m_prm3)
Minimizing 1162      Time: 0:00:00 ( 0.53 ms/it)
false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt          0.04   1.0     .      .      .      .     .
 rcvr         -0.41  -0.47   1.0     .      .      .     .
 trng          0.33   0.34  -0.21   1.0     .      .     .
 drpt & trgt  -0.37   0.25  -0.08  -0.07   1.0     .     .
 rcvr & trgt   0.61   0.04  -0.31   0.13  -0.86   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.3499, 0.5979, 0.7408, 0.8584, 0.9488, 0.9943, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.52   0.03  0.0   -0.01   0.61  -0.57   0.16
 drpt         -0.12   0.62  0.0   -0.06  -0.61  -0.47   0.12
 rcvr          0.33  -0.47  0.0    0.52  -0.19  -0.58  -0.18
 trng         -0.27   0.35  0.0    0.83   0.08   0.31  -0.08
 drpt & trgt   0.45   0.47  0.0   -0.1    0.38  -0.15  -0.63
 rcvr & trgt  -0.58  -0.22  0.0   -0.15  -0.26   0.03  -0.72
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + MixedModels.ZeroCorr((1 + a1 | Item))
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
4: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.5710                          
[2]         41  2356.2006   37.3704      13  0.0004
[3]         42  2341.3452   14.8554       1  0.0001
[4]         56  2353.2482  -11.9030      14     NaN
───────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0760533 0.2757776
drpt 0.0193959 0.1392690 +0.04
rcvr 0.0189445 0.1376392 -0.41 -0.47
trng 0.0350939 0.1873337 +0.33 +0.34 -0.21
drpt & trgt 0.0091654 0.0957362 -0.37 +0.25 -0.08 -0.07
rcvr & trgt 0.0059703 0.0772677 +0.61 +0.04 -0.31 +0.13 -0.86
trgt 0.0002813 0.0167707 . . . . . .
Item (Intercept) 0.0006679 0.0258444
a1 0.0005355 0.0231408 .
Residual 0.0763232 0.2762666

Perhaps not. We stay with LMM m_prm2

4.7 Compare models with goodness-of-fit statistics.

Code
let mods = [m_zcp2, m_prm2, m_prm1, m_cpx];
 DataFrame(;
    model=[:m_zcp2, :m_prm2, :m_prm1, :m_cpx],
    pars=dof.(mods),
    geomdof=round.(Int, (sum  leverage).(mods)),
    AIC=round.(Int, aic.(mods)),
    AICc=round.(Int, aicc.(mods)),
    BIC=round.(Int, bic.(mods)),
  )
end

4 rows × 6 columns

modelparsgeomdofAICAICcBIC
SymbolInt64Int64Int64Int64Int64
1m_zcp228292245024502639
2m_prm242287242524262709
3m_prm148282243024312754
4m_cpx56259246524662843

LMM m_prm2 is a defensible solution according to \(\Delta\) AIC, \(\Delta\) BIC suggests we should not bother with CPs.

Code
coeftable(m_prm2)
Coef. Std. Error z Pr(>
(Intercept) 7.68255 0.0388114 197.95 <1e-99
trgt -0.00948337 0.00442408 -2.14 0.0321
trng -0.0841462 0.0279289 -3.01 0.0026
drpt 0.042879 0.022888 1.87 0.0610
rcvr -0.0370272 0.0226406 -1.64 0.1020
a1 -0.0922378 0.0339009 -2.72 0.0065
trgt & trng -0.00176196 0.00963308 -0.18 0.8549
trgt & drpt -0.0199346 0.0181564 -1.10 0.2722
trgt & rcvr -0.00411183 0.0162151 -0.25 0.7998
trgt & a1 -0.00874189 0.00389127 -2.25 0.0247
trng & a1 0.0414062 0.0244065 1.70 0.0898
drpt & a1 0.0287386 0.0198962 1.44 0.1486
rcvr & a1 -0.0208128 0.0197094 -1.06 0.2910
trgt & trng & a1 -0.00234866 0.00840927 -0.28 0.7800
trgt & drpt & a1 0.017156 0.0158645 1.08 0.2795
trgt & rcvr & a1 0.00435148 0.0141861 0.31 0.7590

4.8 Pruning fixed-effects

There is an option of pruning some higher-order interactions.

4.9 Diagnostic plots

Various visualizations are used to check whether or not data are defensibly modeled with an LMM. They may lead to removal of outliers, transformations of the dependent variable, and deliver valuable heuristic information to be followed up with exploratory post-hoc analyses or ideally replication of new insights gained this way. In practice, it appears that only severe violations will stop people from reporting a model.

4.9.1 Residuals over fitted

Code
scatter(fitted(m_prm2), residuals(m_prm2))

Looks like we missed some fast response times.

4.9.2 Q-Q plot

Code
qqnorm(m_prm2; qqline=:none)

Hm?

4.9.3 Residual distributions: observed vs. theoretical

Curves for residulas based on observed and theoretical values should correspond.

Code
let
  n = nrow(dat)
  dat_rz = (;
    value=vcat(residuals(m_prm2) ./ std(residuals(m_prm2)), randn(n)),
    curve=repeat(["residual", "normal"]; inner=n),
  )
  draw(
    data(dat_rz) *
    mapping(:value; color=:curve) *
    density(; bandwidth=0.1);
  )
end

Figure 1: Kernel density plot of the standardized residuals for model m1 versus a standard normal

They are a bit too narrow.

4.10 Conditional means of random effects

4.10.2 Borrowing-strength plots

Shrinkage refers to the adjustment of subject-level or item-level predictions by taking population estimates into account. The further a subject’s/item’s estimate is from the fixed effect or the more variable or less reliable the subject’s/item’s estimate, the more the prediction will be shrunk towards the population estimate. Alternative terms for shrinkage are “borrowing strength” (Tukey) and regularization. My favorite is actually Tukey’s because indeed we borrow strength from the population estimates to make predictions for individual subjects’ effects. The goal of this section to illustrate the results of borrowing strength.

Subject-related conditional means of random effects revealed information about individual differences beyond fixed effects. Would these results also be visible in unconditional means, that is when we compute GM and experimental effects within subjects (i.e., as fixed effects) without borrowing strength from the population estimates?

In the following plots, effect estimates based on alone on each subject’s data (i.e., no pooling of data, no borrowing of strength) are plotted in pink and the subjects’ conditional means shown in the caterpillar plots are plotted in blue. The arrows indicate how much a subject’s prediction is changed by borrowing strength from knowledge of the population estimates.

Code
shrinkageplot!(Figure(; resolution=(1000, 1200)), m_prm2)

Figure 3: Shrinkage plots of the subject random effects in model m_prm2

There is an outlier subject

Code
cm = raneftables(m_prm2);
sort(DataFrame(cm.Subj), ["(Intercept)", :trng, "drpt & trgt"])
sort(DataFrame(cm.Subj), :trng)

52 rows × 8 columns (omitted printing of 1 columns)

Subj(Intercept)drptrcvrtrngdrpt & trgtrcvr & trgt
StringFloat64Float64Float64Float64Float64Float64
1p10-0.665239-0.08007830.176854-0.5029660.3251-0.281851
2p14-0.5171290.149376-0.104141-0.3804180.0865924-0.0397509
3p760.0792905-0.1004870.10478-0.309599-0.02844710.0204237
4p02-0.622029-0.1970880.260561-0.2516220.0703971-0.126565
5p360.421409-0.03633610.0871325-0.231898-0.1689330.145207
6p390.291872-0.1307340.143494-0.200954-0.01920910.00241466
7p74-0.329558-0.06096160.0915499-0.1793660.0670593-0.0858263
8p670.03832890.00117737-0.0194406-0.155040.0331733-0.0118186
9p06-0.2377130.1392930.0263688-0.117255-0.06677320.0402579
10p540.313212-0.1618670.0464652-0.116672-0.03678070.0410875
11p23-0.0288518-0.07110390.0526805-0.103828-0.1269740.0803409
12p42-0.09737980.288202-0.16646-0.1037430.121444-0.0338319
13p590.4216840.0211085-0.255639-0.102933-0.01229040.0649662
14p45-0.363704-0.2957720.169706-0.09756170.0434683-0.0933397
15p30-0.32999-0.150138-0.0896487-0.0878293-0.0134365-0.00936306
16p83-0.23045-0.0533087-0.160339-0.07021770.0131955-0.0083732
17p84-0.389921-0.2029870.444955-0.0490962-0.1994010.0206631
18p720.176323-0.0684158-0.0524196-0.0387423-0.1194730.117996
19p12-0.148151-0.04484630.0523027-0.0284075-0.06716660.0361057
20p040.138257-0.04476610.0390786-0.027842-0.00322990.00440696
21p770.134281-0.09649470.0522384-0.021795-0.08008060.0499684
22p690.104395-0.0390978-0.0831703-0.01444830.02269670.0349963
23p79-0.0769208-0.01931770.0301171-0.0144253-0.03216780.0103847
24p55-0.198413-0.0638211-0.0784104-0.00928064-0.06245240.0241703
25p270.0253274-0.00351493-0.00407116-0.00610536-0.09761570.0736163
26p280.150447-0.00600531-0.1515170.007017520.0683197-0.0247804
27p250.2676370.137046-0.04742720.00967245-0.0356190.0451562
28p640.2683350.002380780.03385870.0132533-0.0188090.024563
29p570.00282995-0.00253166-0.01534070.01559630.0551781-0.0361016
30p53-0.07836820.02574130.02033310.01802850.0430858-0.0321112

It appears to be Subj “p10”.

5 Appendix

Code
versioninfo()
Julia Version 1.8.1
Commit afb6c60d69a (2022-09-06 15:09 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 8 on 8 virtual cores